function [t, A_bas, A_lo, A_corr, c_bas, c_lo, c_corr] = lo_verify(k0, Nmodes, c0, figs)

if nargin < 4
    close all
    figs = 1;
end

if nargin < 2
    Nmodes = 2;
end

addpath ../
p = params;

t1 = logspace(-6, 0, 200);
dt = max(diff(t1))
t2 = t1(end)+dt:dt:p.Tmax;
t = [0 t1 t2];

t = [0 logspace(-4, log10(p.Tmax), 200)];
% t = linspace(0, p.Tmax, 200);


tau = p.delta * t;
N = length(tau)

dt = 1e-7;
if (dt > min(diff(tau)))
    dt = min(diff(tau))
end

h = -(1 + lambertw(-p.beta / (p.beta - 1) * exp(((-p.beta + p.delta * t) / (p.beta - 1))))) * (p.beta - 1);

k = h * k0;
M = (h - 1 + p.beta) / p.beta * p.Ma_prime;

% compute eigenvalues, eigenfunctions, adjoints
lambda = zeros(Nmodes, N);
v = zeros(p.N, Nmodes, N);
adj = v;
for i = 1:N
    
    
    [ev, tmp, L] = lo_comp_eigs(k(i), M(i));
    
    [ef, ev] = eigs(L, Nmodes, 'SM');
    [tmp, P] = sort(diag(ev), 'descend');
    lambda(:,i) = tmp / h(i)^2;
    v(:,:,i) = ef(:,P);
    
    La = lo_comp_adj(k(i), M(i));
    
    [ef, ev] = eigs(La, Nmodes, 'SM');
    [tmp, P] = sort(diag(ev), 'descend');
    lambda_a(:,i) = tmp;
    adj(:,:,i) = ef(:,P);
    
    
    % normalize eigenfunctions
    for j = 1:Nmodes
        v(:,j,i) = v(:,j,i) / sqrt(ip(v(:,j,i), v(:,j,i), 1, p)) / sqrt(h(i));
        adj(:,j,i) = adj(:,j,i) / ip(adj(:,j,i), v(:,j,i), 1, p) / h(i);
    end
    
    
end

% go through and correct sign flips in the eigenfunctions
for i = 2:N
    for j = 1:Nmodes
        if (v(1,j,i) * v(1,j,i-1) < 0)
            v(:,j,i) = -v(:,j,i);
            adj(:,j,i) = -adj(:,j,i);
        end
        
    end
end



% interpolate eigenvectors for easy differentiation
if Nmodes > 1
    v_i = @(t_i) interp3(1:Nmodes, 1:p.N, tau, v, (1:Nmodes), (1:p.N), t_i);
else
    vv = reshape(v, p.N, N);
    v_i = @(t_i) interp2(tau, 1:p.N, vv, t_i, 1:p.N);
end


% plot the eigenvalues
if figs
    figure
    subplot(2,1,1);
    plot(tau, lambda);
    xlabel('\tau');
    ylabel('Eigenvalues');
    subplot(2,1,2);
    plot(tau, lambda_a);
    xlabel('tau');
    ylabel('Adjoint eigenvalues');
end



% compute numerical solution
z = linspace(0, 1, p.N)';
c0 = 16 * z.^2 .* (1 - z).^2;
num = stab(k0, p, c0);

% extract coeffs
num_i = deval(num, tau / p.delta);
c_num = zeros(Nmodes, N);
for i = 1:N
    for j = 1:Nmodes
        c_num(j,i) = ip(reshape(adj(:,j,i), p.N, 1), num_i(:,i), 1, p);
    end
end

% compute the c_i
c_i = zeros(Nmodes, 1);
for i = 1:Nmodes
    c_i(i) = ip(reshape(adj(:,i,1), p.N, 1), c0, 1, p);
end
fprintf('c_i:\n');
disp(c_i);


% loop again to compute gamma
gamma = zeros(Nmodes, Nmodes, N);
dv = zeros(p.N, Nmodes);


for i = 1:N
    
    
    h_2 = interp1(tau, h, tau(i) + dt);
    v_2 = v_i(tau(i) + dt);
    h_1 = interp1(tau, h, tau(i) - dt);
    v_1 = v_i(tau(i) - dt);
    
    for k = 1:Nmodes
        
        if (i == 1)
            tmp_2 = v_2(:,k);
            dv(:,k) = (tmp_2 - v(:,k,i)) / dt;
        elseif (i == N)
            tmp_1 = v_1(:,k);
            dv(:,k) = (v(:,k,i) - tmp_1) / dt;
        else
            tmp_2 = v_2(:,k);
            tmp_1 = v_1(:,k);
            dv(:,k) = (tmp_2 - tmp_1) / 2 / dt;
        end
        
        if (max(isnan(dv(:,k))) == 1)
            fprintf('NaN detected\n');
            fprintf('i = %d, h_1 = %.4e, h_2 = %.4e, v_1 = %d, v_2 = %d\n', i, h_1, h_2, max(max(isnan(v_1))), max(max(isnan(v_2))));
            return
        end
    end
    
    
    
    for j = 1:Nmodes
        for k = 1:Nmodes
            gamma(j,k,i) = -ip(reshape(adj(:,j,i), p.N, 1), dv(:,k), 1, p);
        end
    end
    
    
end




% compute beta
tmp = zeros(N, Nmodes);
for i = 1:Nmodes
    tmp(:,i) = reshape(gamma(i,i,:), N, 1);
end
beta = cumtrapz(tau', lambda' / p.delta + tmp)';



% compute leading order coeffs
c_lo = zeros(Nmodes, N);
c_bas = c_lo;
for i = 1:Nmodes 
    c_lo(i,:) = c_i(i) * exp(beta(i,:));
    c_bas(i,:) = c_i(i) * exp(cumtrapz(tau', lambda(i,:)' / p.delta)');
end

fprintf('c_lo(t = 0):\n');
disp(c_lo(:,1));
fprintf('c_num(t = 0):\n');
disp(c_num(:,1));

% do the corrections
c_corr = zeros(Nmodes, N);
for j = 1:Nmodes
    for l = 1:Nmodes
        if (l ~= j)
            int = cumtrapz(tau, reshape(gamma(j, l,:) .* gamma(l, j,:), 1, N) ./ (lambda(j,:) - lambda(l,:)));
            
            tmp = -c_i(l) * reshape(gamma(j, l, :), 1, N) ./ (lambda(j, :) - lambda(l, :)) .* exp(beta(l,:)) + ...
                (c_i(j) * int + gamma(j, l, 1) * c_i(l)  / (lambda(j,1) - lambda(l,1))) .* exp(beta(j,:));
                        
            c_corr(j,:) = c_corr(j,:) + tmp;
        end
    end
end

% plot the amplitudes
A_num = max(abs(num.y));
A_lo = zeros(1,N);
A_corr = A_lo;
A_bas = A_lo;
for i = 1:N
    A_bas(i) = max(abs( reshape(v(:,:,i), p.N, Nmodes) * c_bas(:,i)));
    A_lo(i) = max(abs( reshape(v(:,:,i), p.N, Nmodes) * c_lo(:,i)  ));
    A_corr(i) = max(abs( reshape(v(:,:,i), p.N, Nmodes) * (c_lo(:,i) + p.delta * c_corr(:,i)) ));
end

if figs
    figure
    plot(num.x * p.delta, A_num / norm(c0, 'inf'), 'k', tau, A_lo / norm(c0, 'inf'), 'b-.', tau, A_corr / norm(c0, 'inf'), 'r--');
    xlabel('\tau');
    ylabel('Amplitudes');
    l = legend('numerical','leading order','leading plus correction', 'location', 'best');
    figure
    loglog(num.x * p.delta, A_num / norm(c0, 'inf'), 'k', tau, A_lo / norm(c0, 'inf'), 'b-.', tau, A_corr / norm(c0, 'inf'), 'r--');
    xlabel('\tau');
    ylabel('Amplitudes');
    l = legend('numerical','leading order','leading plus correction', 'location', 'best');
end

% rescale amps
A_bas = A_bas / norm(c0, 'inf');
A_lo = A_lo / norm(c0, 'inf');
A_corr = A_corr / norm(c0, 'inf');


% plot the coefficients
if figs
    figure
    for j = 1:Nmodes
        subplot(ceil(Nmodes/2),2,j);
        semilogx(tau, c_num(j,:), 'k', tau, c_lo(j,:), 'b-.', tau, c_lo(j,:) + p.delta * c_corr(j,:), 'r--');
        xlabel('\tau');
        ylabel(strcat('c_',num2str(j)));
        tmp = get(gca, 'xlim');
        xlim([tmp(1), p.Tmax * p.delta]);
    end
end